%CES Version

%This file computes the gains in GE from all countries implemeting the
%optimal industrial policies at the same time (Table 3, column 3 in BCDR).


clear all;


K=32; %number of sectors
N=62; %includes ROW
K_man=15;    %number of manufacturing sectors

%%%Load data
%%%Data is organized as follows: first column origin country codes, second column
%%%sector codes, third column sectoral value added, fourth column aggregate origin trade balance,
%%%fifth column origin-sector expenditurs, 6th column starts the bilateral
%%%trade flows



%Note: data is sorted by sector, then origin.

A=importdata('simulation_ge_data.csv');
B=importdata('sim_params_no_int.csv'); %these are the estimated elasticities

%Data vectors
X_ink = A(:,6:5+N); %Matrix of trade flows, (NxK)xN
D_i=A(1:N,4); %vector of trade balances, Nx1
V_ik=A(:,3); %vector of sectoral value added, (NxK)x1
expenditures_ik=A(:,5); %vector of expenditures, (NxK)x1
countryid=A(1:N,1); %vector of country ids, Nx1

%Parameter vectors
theta_nm=4.32; %non-manufacturing sector TE set to median of manufacturing
h_theta = theta_nm*[ones(2,1);zeros(K_man,1);ones(K-2-K_man,1)]+[zeros(2,1);B(:,2);zeros(K-2-K_man,1)]; %trade elasticities
h_gamma = [zeros(2,1);B(:,1);zeros(K-2-K_man,1)]; %Heterogeneous scale elasticities, non-manufacturing set to zero
gammaXtheta=h_theta.*h_gamma;
rho=.87; %Upper tier elasticity of substitution

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reshaping data vectors

%Reshape trade matrix to be NxNxK (3 dimensional array)
%rows are exporters, columns are importers, 3rd dimension is sectors
X=zeros(N,N,K);
for k=1:K
    X(:,:,k) = X_ink(1+(k-1)*N:k*N,:);
end
X_ink=X;

%Other data matrices
V_ik = reshape(V_ik, [N 1 K]);
D_i =reshape(D_i, [N 1 1]);
expenditures_ik = reshape(expenditures_ik,[N,1,K]);
expenditures_nk=permute(expenditures_ik,[2 1 3]);

%Create input matrices
expenditures_nnk=repmat(expenditures_nk,[N 1 1]);
lambda_ink = X_ink./expenditures_nnk; %trade shares

expenditures_n=sum(expenditures_nk,3);
beta_nk=expenditures_nk./repmat(expenditures_n,[1 1 K]);

V=sum(sum(V_ik));
V_ik=V_ik/V; %normalizing by world value added
D_i = D_i/V; %normalizing by world value added


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gains from Industrial Policy

%Defining output vectors
gfip_hte_0=zeros(N,1); %Het TE, Het SE (service 0 se)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gains from trade and industrial policy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Computing gains from all countries implementing policies at the same
%time

%%%%%%%%%%%%
%Het TE, Het SE (service 0 SE)
theta=h_theta;
gamma = h_gamma;
gamma=reshape(gamma,[1 1 K]);
gamma = repmat(gamma,[N 1 1]);
theta=reshape(theta, [1 1 K]);

%Now industrial policy only
s_k=gamma;
s_ink=repmat(s_k,[1,N,1]);
t_ink=zeros(N,N,K);
[model_output]=gfop_ces_ge_solver(V_ik,D_i,lambda_ink,beta_nk,gamma,theta,rho,t_ink,s_ink,K,N);
gfip_hte_0=model_output.welfare_i;

%Saving the gains
gains=[countryid gfip_hte_0];
csvwrite('gains_hte_ces_ge.csv',gains);

